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GAS  CHANNEL  FORMATION  STUDIES  INCLUDING 
REAL-AIR  THERMODYNAMIC  PROPERTIES 


1.  Introduction 

Hot  reduced-density  channels  can  be  formed  in  gases  by  methods  such 

12  3 
as  exploding  wires  ,  lasers  ,  pulsed  electrical  discharges  ,  or  charged 

4 

particle  beams  .  In  many  cases,  the  heating  occurs  very  rapidly  com¬ 
pared  to  the  hydrodynamic  response  time  of  the  gas  (the  sound  transit 
time  over  the  channel  radius),  and  can  thus  be  treated  as  instantaneous. 
Furthermore,  cooling  due  to  heat  conduction  and  radiation  is  slow,  for 
many  applications,  and  can  be  neglected  in  treating  the  hydrodynamic 
response.  Under  these  conditions,  channel  formation  is  a  one¬ 
dimensional  hydrodynamic  process  which  departs  from  adiabaticity  only 
because  of  the  occurrence  of  shock  waves.  Analytic  solutions  for  this 
process,  in  the  form  of  self-similar  blast  waves5,  are  appropriate  at 
very  high  initial  overpressure,  but  in  the  range  of  primary  interest 
to  us,  with  overpressure  less  than  a  hundred  times  ambient,  it  is 
necessary  to  run  one-dimensional  fluid  codes  to  calculate  the  time 

evolution  of  the  overheated  channel. 

£ 

We  have  previously  reported  fluid  code  studies  of  this  process 
for  the  case  in  which  the  heating  source  has  a  very  smoothly  rounded 
radial  profile,  e.g.  the  Bennett  profile  appropriate  for  a  self-pinched 
electron  beam.  In  this  case,  we  found  that  the  channel  formation 
process  went  through  the  following  stages:  (i)  Initially,  the  gas  is 
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at  uniform  density,  but  with  an  overpressure  profile  corresponding  to 
the  source  heating  profile,  (ii)  A  rarefaction  wave  reduces  the  gas 
density  within  the  heating  source  radius,  (iii)  An  outgoing  shock  wave 
forms  near  the  point  of  steepest  slope  of  the  initial  overpressure 
profile,  and  leads  to  expansion  of  the  gas  over  a  channel  broader  than 
the  source  radius.  (iv)  When  the  shock  detaches  from  the  channel,  the 
entire  channel  is  left  at  a  pressure  of  about  85-907.  of  ambient. 

(v)  The  channel  slowly  returns  to  ambient  pressure,  with  a  density 
rise  of  about  10%  and  no  further  shock  formation.  (vi)  The  final 
density  profile  is  determined  strictly  by  adiabatic  expansion  in  the 
central  channel ,  but  the  density  is  further  depressed  in  the  radial 
wings  by  the  entropy  increase  attendant  upon  shock  heating.  Under 
these  conditions,  we  gave  an  analytic  formulation  of  the  density 
profile  at  late  times,  when  pressure  balance  has  been  restored. 6  This 
formulation  depended  on  the  "polytropic"  assumption  that  the  ratio  y  of 
the  specific  heat  c^  (at  constant  pressure)  and  cv  (at  constant  volume) 
is  a  constant,  which  we  usually  took  to  be  y  =  cp/cv  *  1.4,  charac¬ 
teristic  of  a  diatomic  gas  at  moderate  temperature.  For  simplicity 
and  comparability,  we  used  the  same  assumption  in  our  fluid  code  runs. 

In  a  real  diatomic  gas,  y  depends  on  both  temperature  T  and  density 
and  decreases  from  'v*  1.4  at  room  temperature  to  less  than  1.2  at 
T  >  3000°K.  The  purpose  of  the  present  rate  is  to  include  this 
variation  in  the  analytic  formulation,  which  becomes  possible  because 
a  generalized  form  of  y  can  be  expressed  as  a  function  of  specific 
energy  e,  with  practically  no  further  dependence  on  p .  We  also  show 
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the  effect  of  including  real  gas  thermodynamics  in  the  fluid  code 
runs. 

In  another  previous  report7,  we  have  shown  that  the  hydrodynamic 

evolution  becomes  somewhat  more  complex,  with  ringing  of  the  channel 

and  sometimes  multiple  shocks  occurring,  if  the  radial  profile  of  the 

heating  source  is  squarer  in  shape,  or  has  discontinuities.  These 

complications  are  treated  only  in  passing  in  the  present  note. 

In  Sec.  2,  we  give  a  convenient  summary  and  formulation  of  the 

real  gas  thermodynamic  properties  required.  In  Sec.  3,  we  include 

these  properties  in  a  semi-analytic  model  of  the  density  profile  at 

late  time.  In  Sec.  4,  we  compare  the  results  of  the  model  with 

numerical  studies  performed  on  the  one-dimensional  (cylindrical) 

6—9 

Eulerian  flux-corrected  transport  fluid  code  FAST1D. 
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2.  Real  Gas  Thermodynamics 

A  particularly  convenient  way  of  specifying  the  thermodynamic 
properties  of  a  real  gas  is  by  means  of  a  constitutive  equation. 

P/p  -  ef (e,p) ,  (1) 

supplemented  by  an  ideal  gas  equation  of  state, 

P/P  -  (kB/M)[l  +  a(e,p)]T,  (2) 

where  P  is  pressure,  p  is  density,  e  is  specific  energy  (thermal  plus 
internal,  but  excluding  any  macroscopic  kinetic  energy) ,  T  is  tempera¬ 
ture,  k  is  Boltzmann's  constant,  M  is  the  molecular  mass,  and  f  and  a 
B 

are  functions  that  characterize  the  gas  species.  The  quantity  a 
measures  the  degree  of  dissociation  and  ionization;  it  is  unity  for 
the  undissociated,  unionized  state,  and  is  equal  in  general  to  the  mean 
number  of  dissociated  atoms,  ions  and  electrons  per  molecule.  The 
quantity  f  is  a  generalization  of  the  adiabatic  index  y ,  defined  by 

y  =  c  /c  ,  (3) 

p  v 

for  the  case  in  which  the  specific  heats  c^  (at  constant  pressure)  and 
cv  (at  constant  volume)  are  constant.  In  this  case, 

f(e,P)  =  Y  -  1,  (4) 

and  the  adiabatic  relations 

P/pY  *  const  ,  (5) 

e/pY  1  *  const  (6) 

can  be  used  to  characterize  adiabatic  flows. 
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In  a  diatomic  gas  such  as  air,  f(e,p)  ranges  from  2/3  at  cryogenic 
temperatures  where  e  consists  entirely  of  translational  thermal  energy, 
to  2/5  at  room  temperature,  where  rotational  states  are  fully  excited, 
to  2/7  at  T  ^  3000°K,  where  vibrational  states  are  fully  excited.  Up 
to  this  point,  f(e,p)  is  quite  independent  of  density.  At  still  higher 
temperatures,  f(e,p)  decreases  further  as  electronic  excitation,  disso¬ 
ciation,  and  ionization  begin  to  contribute  to  e.  Since  the  latter  two 
processes  change  the  number  of  particles,  they  lead  to  a  p  dependence 
of  f.  However  the  particular  form  (1)  was  chosen  because  this  p 
dependence  remains  minimal,  even  in  this  regime.  This  point  is 
illustrated  in  Fig.  1,  where  f(e,p),  taken  from  the  tabulated  values 
of  Ref.  10,  is  plotted  for  air  at  a  variety  of  densities,  at  tempera¬ 
tures  up  to  ^  10,000ok,  and  it  is  seen  that  varying p  over  a  wide  range 
changes  f(e,p)  no  more  than  10%.  Consequently,  we  can  regard  f  as  a 
function  of  e  alone, 

f(e,P)  =  f(t),  (7) 

which  leads  to  considerable  simplification  in  our  analysis.  If  f  were 
regarded  as  a  function  of  T  and  p ,  or  if  a  relation  were  written 
between  e  and  T,  the  p  dependence  would  be  much  stronger.  A  second 
advantage  of  the  formulation  used  in  Eqs.  (1)  and  (2)  is  that,  for 
most  purposes,  we  do  not  need  to  know  T,  so  that  only  Eq.  (1)  is  used. 
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Density  Profile  at  Late  Times 
The  first  law  of  thermodynamics, 

Tds  =  de  +  Pdp"1  ,  (8) 

can  be  rewritten  in  the  form 


dp  _  de  _  M  ds 
p  ef(e)  k  1+a 

B 


(9) 


where  s  is  the  specific  entropy  and  Eqs.  (1)  and  (2)  have  been  used. 
To  find  the  density  profile  at  late  times,  when  P  =  P  is  uniform  at 

cl 

the  ambient  value,  we  can  integrate  (9)  along  the  trajectory  of  an 
element  of  gas,  to  obtain 


°-l££±  .  exp 


as  (r,00) 
e (r  ,q) 


de  _ 

s  *f  (e' ) 


(10) 


where  an  element  of  gas  at  r  initially  moves  to  r  at  t  *  00 ,  A  is  the 

o  s 

integral  of  the  second  term  in  (9) ,  and  we  have  used  the  fact  that 
p  =  p  is  uniform  at  t  =  0.  In  Eq.  (10),  the  value  of  e(r,00)  is  deter- 
mined  by  applying  Eq.  (1)  : 


Pa/p(r,«)  =  e  (r ,»)  f  [e  (r,08)]  .  (11) 

Usually  the  energy  deposition  profile  e (r  ,0)  is  given;  it  is  related 

o 

to  the  initial  pressure  profile  by 


P(r  ,0)/p  *  e(r  ,0)f[e(r  ,0)1. 

o  a  o  *■  o  -* 


(12) 
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Finally,  r  is  related  to  r  by  integrating  the  relation 

o 

P (r  ,0)  2 

o  _  dr  _ 

P(r,°°)  dr  2  ' 

o 


to  obtain 


dr  2r  '  p(r' ,”) . 


(13) 


Equations  (10)  -  (13),  along  with  the  tabulated  function  f(e)  from 

Fig.  1,  determine  the  density  profile  at  t  =  if  Ag  can  be  specified. 

An  entropy  change  occurs  only  when  the  outgoing  shock  wave  passes 

over  a  particular  element  of  gas.  This  does  not  occur  at  all  in  the 

£ 

central  region  of  the  channel  ,  and  only  once  in  the  outer  range  of  r. 
Thus 


A 

s 


u 


_d§_ 


1  +  a 


where  the  integral  is  taken  over  the  shock  front,  for  conditions 
prevailing  when  the  shock  passes  over  a  particular  element  of  gas. 

Shock  heating  typically  results  in  a  fairly  small  change  in  the 
final  density  profile,  so  a  lowest  order  approximation  is  given  by 
simply  neglecting  As»  In  Reference  6,  we  went  one  step  further  by 
giving  a  heuristic  formula  for  A^r^)  for  the  case  of  a  Bennett  heating 
profile,  based  on  empirical  results  of  fluid  code  runs  with  y  -  7/5,  as 
well  as  the  Rankine-Hugoniot  relations  and  certain  other  general  formal 


requirements. 


We  see  little  hope  of  writing  down  any  improved  closed  form  for  that 
applies  to  the  more  complex  situation  considered  here,  where  y  is  no 
longer  constant.  However  the  effect  of  the  shock  is  merely  a  small  to 
moderate  correction  to  the  density  profile  at  late  times,  and  further¬ 
more  the  shock  effect  is  greatest  in  the  radial  wings,  where  the  gas  is 
not  too  hot  and  the  assumption  that  y  =  7/5  is  closer  to  reality.  Thus 
we  shall  simply  continue  to  use  Eqs.  (15)  and  (16)  in  Eq.  (10)  to 
specify  the  density  profile.  As  we  shall  see,  the  results  of  this 
procedure  are  quite  acceptable  when  compared  to  fluid  code  runs. 


4. 


Results 


It  is  most  important  to  use  Eq.  (1)  with  the  correct  dependence 
f(e)  to  calculate  the  initial  overpressure  profile  p(r,0)  resulting 
from  a  given  profile  e(r,0)  of  energy  deposited  by  the  source.  As 
seen  in  Fig.  1,  the  pressure  can  be  overestimated  by  a  factor  of  two 
or  more  if  the  assumption  y  =  7/5,  i.e.  f(e)  =  0.4,  is  used  to  calcu¬ 
late  P(r,0).  In  more  familiar  terms,  the  specific  heat  increases 
significantly  with  temperature  [but  it  is  more  convenient  to  work  with 
the  function  f (e)  defined  in  Eq.  (1) ,  rather  than  with  the  specific 
heat.]  The  shape  as  well  as  the  amplitude  of  the  pressure  profile  is 
altered;  in  general,  P(r,0)  will  be  broader  and  more  square-shaped 
than  e(r,0),  as  seen  in  Figs.  2  and  3. 

In  Figs.  2-9,  the  gas  is  taken  to  be  air  at  standard  temperature 
at  pressure  before  heating  by  the  external  source.  The  energy  deposi¬ 
tion  profile  of  the  external  source  is  taken  to  be  Bennett,  i.e.  e(r, 
t  =  0)  -  e  has  a  Bennett  shape  (where  z  is  the  ambient  specific 

3  3. 


energy) , 


C(r't  =  0)  =  £a  *  d*l°'33/W 


Dimensionless  units  for  r  and  t  are  used  in  the  figures,  r  =  r/a,  and 


t  =  c  t/a,  where  c  is  the  ambient  sound  speed.  We  also  use  dimen- 
sa  sa 

sionless  units  for  P  and  p,  defined  by  ?  =  P/P  and  p  =  p/p  ,  where 

a  a 

P  and  p  are  the  ambient  values.  The  overpressure  factor  P  ,  defined 

3  3  O 

by 


Pq  =  P(?  =  0,  t  =  0)  -  1, 


is  frequently  used  to  characterize  the  strength  of  the  initial  energy 


deposition.  For  the  polytropic  model  which  we  studied  previously, 
this  choice  of  units  reduced  the  problem  to  completely  dimensionless 
form;  this  is  not  true  in  the  real  gas,  however,  since  f(e)  depends 
on  the  specific  energy  density  e  and  cannot  be  scaled  out. 

m  order  to  examine  the  effect  of  real-gas  thermodynamics  on  the 

hydrodynamic  flow  subsequent  to  energy  deposition,  we  shall  compare 

cases  in  which  e(r,0)  has  a  Bennett  profile,  and  (i)  the  assumption 

y  -  7/5  is  used  to  calculate  P(r,0)  and  also  to  treat  the  subsequent 

expansion,  and  the  amplitude  of  e(r,0)  is  such  that  Pq  has  same  given 

value;  (ii)  the  accurate  dependence  of  f(e)  is  used  to  calculate 

P(r,0)  and  also  to  treat  the  subsequent  expansion,  and  the  amplitude 

of  e(r,0)  is  chosen  such  that  P  has  the  same  as  in  case  (i) .  Case 

o 

(i)  will  be  referred  to  as  the  polytropic  case,  (ii)  as  the  real  gas 
case. 

Figure  4  shows  snapshots  of  the  density  profile  p (r)  at  a  series 
of  times,  for  the  real  gas  case  with  Pq  =  11.5.  The  phenomenology  of 
channe 1  formation  —  initial  rarefaction  wave,  deepening  and  broadening 
of  the  channel,  shock  formation  and  propagation,  density  undershoot 
at  late  time ,  and  slow  return  to  pressure  balance  —  is  qualitatively 
similar  to  that  seen  with  the  constant  y  assumption.^  Fig.  5  compares 
the  pressure  profile  at  a  particular  time  for  real  gas  and  polytropic 
cases.  We  see  that  the  shock  is  a  little  stronger  for  the  real  gas 
than  for  the  corresponding  polytropic  case.  To  separate  the  effect 
of  the  broadened  P(r,0)  profile  in  the  real  gas  case  from  the  effect 
of  the  variation  of  f(e)  during  expansion,  we  ran  another  case  in 


which  the  P(r,0)  profile  was  calculated  from  the  exact  thermodynamics, 
but  f(e)  was  set  equal  to  0.4  during  the  subsequent  expansion.  In 
this  case,  shown  in  Fig.  6,  the  shock  was  a  little  weaker  for  the  real 
gas  case.  Thus  both  the  change  in  the  P(r,0)  profile  and  the  varia¬ 
tion  of  f(e)  during  the  subsequent  expansion  have  some  influence  on 
the  phenomenology  (and  in  opposite  directions) ,  but  the  differences 

from  the  ideal-gas  case  are  small,  for  a  given  value  of  P  . 

o 

In  Figs.  7  and  8,  we  compare  the  time  dependence  of  the  on-axis 

density  p(0,t)  for  the  real  gas  and  polytropic  cases.  We  see  that  in 

the  polytropic  case,  the  p(0,t)  curve  is  smooth,  dips  slightly  below 

the  final  equilibrium  value  and  then  returns  to  equilibrium  slowly 
6 

and  monotonically  ,  but  in  the  real  gas  case,  the  channel  density 

7 

exhibits  weak  oscillations.  We  have  shown  in  previous  work  that  this 
type  of  behavior  is  characteristic  of  pressure  profiles  that  are 
squarer-shaped  than  Bennett;  thus  it  would  appear  to  be  due  to  the 
squaring  up  of  the  P(r,0)  profile  that  results  when  the  correct  f(e) 
is  used  to  calculate  P(r,0).  However  the  effect  is  again  small,  for 
practical  purposes. 

Finally,  in  Fig.  9,  we  compare  the  results  of  the  model  developed 
in  Sec.  3,  Eqs.  (10)  -  (15)  with  p (r)  profiles  observed  at  late  times 
in  a  code  run.  We  find  that  the  model  gives  a  good  representation  of 
the  density  profile  at  late  times,  for  all  cases  that  have  been  run. 


5.  Conclusions 

It  is  important  to  use  accurate  values  of  the  specific  heat,  or  of 
the  function  f(e),  in  calculating  the  initial  pressure  profile  asso¬ 
ciated  with  a  given  energy  deposition  source.  Changes  in  the  function 
f(e)  in  air  at  a  few  thousand  degrees  reduce  the  initial  overpressure 
by  about  a  factor  of  two,  and  also  make  the  initial  pressure  profile 
considerably  squarer-shaped  than  the  energy  deposition  profile. 

However  the  effect  of  using  accurate  values  of  f(e)  rather  than 
a  constant  value  y  =  7/5  during  the  hydrodynamic  expansion  is  small, 
with  the  most  noticeable  effect  being  the  slight  ringing  of  the  on-axis 
density  due  to  the  squarer  shape  of  the  P(r,0)  profile.  For  most 
practical  purposes,  it  is  adequate  to  use  y  =  7/5  during  the  expansion, 
provided  P(r,0)  is  calculated  accurately.  If  greater  accuracy  is 
needed,  it  is  convenient  to  formulate  the  thermodynamic  properties  of 
the  gas  in  terms  of  the  function  f{e),  rather  than  through  specific 
heats  or  other  temperature-dependent  quantities.  The  semi-analytic 
model  of  Sec.  3  gives  an  accurate  representation  of  the  density  profile 
of  the  heated  channel  at  late  times,  when  pressure  balance  has  been 
restored. 
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Fig.  1.  Thermodynamic  properties  of  air.  The  function  f(c)  defined 


in  Eq.  (1)  is  plotted  against  specific  internal  energy  e,  at 

various  densities:  dots  correspond  to  p  •  7.0  x  10  5  gm/cm3 

-4  3 

(0.054  atm),  squares  to  p  »  3.0  x  10  gm/cm  (0.23  atm), 
crosses  to  o  ■  1.29  x  10  3  gm/cm3  (1  atm). 
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6.  Radial  density  profile  at  "t  =  6.33  for  two  cases  with 

PQ  ■  26.1  and  Bennett  energy  deposition  profile.  In  each 
case,  the  initial  pressure  profile  P(r,0)  is  calculated 
using  the  real  gas  function  f(e).  For  Curve  I,  the 
polytropic  model  with  y  =  7/5  is  used  to  treat  the  subsequent 
expansion;  for  Curve  II,  the  real  gas  properties  are  used. 
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initial  overpressure  factor  P  =  26.1,  in  real  gas. 


Radial  density  profiles  at  late  times  in  real  gas.  Curve  I 
analytic  solution.  Curve  II:  hydrodynamic  code  solution  at 
t  =  9.66.  Curve  III:  hydro  code  solution  at  t  =  18.0.  The 
initial  overpressure  factor  was  P  =11.5. 
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